
#pragma OPENCL EXTENSION all : enable

#define ALIGN     16
#define ROUND(n)  (((size_t)(n)+ALIGN-1)&(~(size_t)(ALIGN-1)))
#define LOAD(i, Y, j, N) { i.x = Y[j]; i.y = Y[j + N]; i.z = Y[j + 2*N]; }

// This kernel is executed with a very high number of threads = total number of vertices * 3
__kernel
void eating(__global float *vbo, float dy, unsigned vertices_per_atom)
{
   int index = get_global_id(0);
   int tpa = 3 * vertices_per_atom; // threads per atom
   // ghost or pacman?
   float enable = (index / tpa) % 2 ? 0.0f : 1.0f;
   // y coordinate?
   float is_y = ((index % 3) == 1) ? 1.0f : 0.0f;
   // vertex belongs to upper part of sphere?
   float is_up = (index % tpa >= tpa/2) ? 1.0f : -1.0f;

   vbo[index] += dy * is_y * is_up * enable;
}

__kernel
void move_vertices(__global float *vbo, __global float *speed, unsigned vertices_per_atom)
{
  int index = get_global_id(0);

  // threads par atome
  unsigned tpa = (3*vertices_per_atom);

  // atomes au total
  unsigned int Ntot = ROUND(get_global_size(0) / tpa);

  // index de cet atome
  unsigned int Nidx = index / tpa;
  
  vbo[index]= vbo[index]+speed[(index % 3)*Ntot + Nidx];
}

__kernel
void update_position(__global float *pos, __global float *speed, unsigned N)
{
  int index = get_global_id(0);
  pos[index] = pos[index] + speed[index];
}

__kernel
void border_collision(__global float *pos, __global float *speed, __constant float *min, __constant float *max, float radius, unsigned N)
{
  int index =  get_global_id(0);
  if ( pos[index] > ( max[index/ROUND(N)] - radius) || pos[index] < ( min[index/ROUND(N)] + radius) ) {
     speed[index] = -speed[index];
  } 

}

void choc_elastique(int start, int indice,__global float *pos, __global float *speed, int r) {
	float3 Ca = {pos[start],pos[start + r],pos[start + (2*r)]};
	float3 Cb = {pos[indice],pos[indice + r],pos[indice + (2*r)]};
	float3 i =  normalize(Cb - Ca);
	
	float3 Va = {speed[start],speed[start + r],speed[start + (2*r)]};
	float3 Vb = {speed[indice],speed[indice + r],speed[indice + (2*r)]};
	
	float3 j = {0,-i.z,i.y};
	j = normalize(j);
	
	//float3 k = {i.y*j.z - i.z*j.y,-(i.x*j.z - i.z*j.x),i.x*j.y - i.y*j.x};
	float3 k = {(i.y*j.z) - (i.z*j.y),-(i.x*j.z),(i.x*j.y)};
	float3 M[3] = { 
    i,
	  {-i.y, i.x+(pow(i.z, 2) / (1+ i.x)), (-i.y*i.z) / (1+i.x)},
	  {-i.z, (-i.y*i.z) / (1+i.x), i.x + (pow(i.y, 2)/(1 + i.x))}
	};
	
	float3 Var =  {dot(M[0],Va),dot(M[1],Va),dot(M[2],Va)};
	float3 Vbr =  {dot(M[0],Vb),dot(M[1],Vb),dot(M[2],Vb)};
	
	float3 Vta = {Vbr.x,Var.y,Var.z};
	float3 Vtb = {Var.x,Vbr.y,Vbr.z};
	
	float3 temp[3] = {
    {
      M[0].x,
      M[1].x,
      M[2].x
    },
    {
      M[0].y,
      M[1].y,
      M[2].y
    },
    {
      M[0].z,
      M[1].z,
      M[2].z
    }
  };
  
	float3 Vaf = {dot(temp[0],Vta),dot(temp[1],Vta),dot(temp[2],Vta)};
	float3 Vbf = {dot(temp[0],Vtb),dot(temp[1],Vtb),dot(temp[2],Vtb)};
	
	speed[start] = Vaf.x;
	speed[start + r] = Vaf.y;
	speed[start + (2*r)] = Vaf.z;
	
	speed[indice] = Vbf.x;
	speed[indice + r] = Vbf.y;
	speed[indice + (2*r)] = Vbf.z;
	
	
} 


// Version 1
/*
__kernel
void atom_collision(__global float *pos, __global float *speed, float radius)
{
  int index = get_global_id(0);
  int nb_thread = get_global_size(0);
  int i = 0;
  int r = ROUND(nb_thread);
  float3 pos_index = {pos[index], pos[index + r], pos[index + (2*r)]};
  
  for (i= index+1;i<nb_thread;i++) {
    float3 pos_tmp = {pos[i], pos[i + r], pos[i + (2*r)]};
    if (distance(pos_index,pos_tmp) <= 2*radius ) {
      choc_elastique(index,i,pos,speed,r);
    }
  }
} 
*/
// Version 2
/*
void atom_collisionv2_aux(__global float *pos, __global float *speed, int start, int nb_atoms, float radius) {
  int i, j;
  int r = ROUND(nb_atoms); 
  
  float3 pos_start = {pos[start], pos[start + r], pos[start + (2*r)]};
  for (i = start + 1; i < nb_atoms; i++) {
    float3 pos_tmp = {pos[i], pos[i + r], pos[i + (2*r)]};
    if (distance(pos_start,pos_tmp) <= 2*radius ) {
      speed[start] = 0;
      speed[i] = 0;
      speed[start + r] = 0;
      speed[i + r] = 0;
      speed[start + (2*r)] = 0;
      speed[i + (2*r)] = 0;
    }
  }
}

__kernel
void atom_collision(__global float *pos, __global float *speed, float radius)
{
  int index = get_global_id(0);
  int nb_atoms = get_global_size(0) * 2;
  atom_collisionv2_aux(pos, speed, index, nb_atoms, radius);
  atom_collisionv2_aux(pos, speed, nb_atoms - index - 1, nb_atoms, radius);
} 

*/

/*
void atom_collisionv3_aux(
    __global float *pos,
    __local float3 *pos_tmp, 
    __global float *speed, float radius, 
    int start_group,
    int start_id
    ) {
  int nb_atoms = get_global_size(0) * 2;
  int i, j;
  int r = ROUND(nb_atoms);
  int lid = get_local_id(0);

  float3 pos_index = { pos[start_id], pos[start_id + r], pos[start_id + (2*r)] };
 
  for (i = start_group; i < get_num_groups(0) * 2; i++) {
    int idx = i * get_local_size(0) + get_local_id(0);
    pos_tmp[lid].x = pos[idx];
    pos_tmp[lid].y = pos[idx + r];
    pos_tmp[lid].z = pos[idx + (2*r)];
    barrier(CLK_LOCAL_MEM_FENCE);

    for (j = 0; j < get_local_size(0); j++) {
      int idx = i * get_local_size(0) + j;
      if (idx != start_id && distance(pos_index, pos_tmp[j]) <= 2*radius) { 
         speed[start_id] = 0;
         speed[idx] = 0;
         speed[start_id + r] = 0;
         speed[idx + r] = 0;
         speed[start_id + (2*r)] = 0;
         speed[idx + (2*r)] = 0;
      }
    }
  }
} 
*/

#define ATOM_COLLISION_WG_SIZE 16
__kernel
void atom_collision(__global float *pos, __global float *speed, float radius)
{
//  Version 3:
//   __local float3 pos_tmp[ATOM_COLLISION_WG_SIZE];
//  int lid = get_local_id(0), gr_id = get_group_id(0);
//  int gl_id = get_global_id(0);
//  int nb_atoms = get_global_size(0) * 2;
//  int nb_groups = get_num_groups(0) * 2;
//  atom_collisionv3_aux(pos, pos_tmp, speed, radius, 
//    gr_id, gl_id);
//  atom_collisionv3_aux(pos, pos_tmp, speed, radius, 
//    nb_groups - 1 - gr_id,
//    nb_atoms - 1 - gl_id
//    );
//
//  

  // Version 3b:
  __local float3 pos_a[ATOM_COLLISION_WG_SIZE];
  __local float3 pos_b[ATOM_COLLISION_WG_SIZE];

  int lid = get_local_id(0);
  int agid = get_global_id(0);
  int ls = get_local_size(0);
  int gs = get_global_size(0);
  int r = ROUND(gs);
  int i, j;

  LOAD(pos_a[lid], pos, agid, r);

  for (i = 0; i < get_num_groups(0); i++) {
    int bgid = i * ls + lid;
    LOAD(pos_b[lid], pos, bgid, r);
    barrier(CLK_LOCAL_MEM_FENCE);

    for (j = 0; j < ls; j++) {
        bgid = i * ls + j;
        if (bgid != agid && distance(pos_a[lid], pos_b[j]) <= 2*radius) {
          choc_elastique(agid,bgid,pos,speed,r);
        }
    }
  }
}
__kernel
void gravity(__global float *pos, __global float *speed, float g)
{
	int index =  get_global_id(0);
	speed[index + ROUND(get_global_size(0))] -= g;
}


void lennard_jones_aux(
    int start, int indice,__global float *pos,
    __global float *speed, int r, int radius
    ) {

  float3 Ca = {pos[start],pos[start + r],pos[start + (2*r)]};
  float3 Cb = {pos[indice],pos[indice + r],pos[indice + (2*r)]};
  float3 i =  normalize(Cb - Ca);
  float d =  distance(Cb, Ca);
  float sigma = 1.25;
  float epsilon = 0.01;
  
  float intensity = 4*epsilon*(12*pow(sigma,12.0f)/pow(d, 13.0f) 
                               - 6*pow(sigma, 6.0f) / pow(d, 7.0f));
  float3 force = {intensity * i.x, intensity * i.y, intensity * i.z};
  
  speed[start] += force.x;
  speed[start + r] += force.y;
  speed[start + (2*r)] += force.z;
  
  speed[indice] -= force.x;
  speed[indice + r] -= force.y;
  speed[indice + (2*r)] -= force.z;
}

__kernel
void lennard_jones(__global float *pos, __global float *speed, float radius)
{
  int index = get_global_id(0);
  int nb_thread = get_global_size(0);
  int i = 0;
  int r = ROUND(nb_thread);
  
  for (i = index + 1; i < nb_thread; i++) {
    lennard_jones_aux(index, i, pos, speed, r, radius);
  }
}




